;************************************************
; NCL Graphics: 
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"

;************************************************
begin

  load "load_setup.ncl"

  vname  = "homhet"
  vunits = "#/L" 
  vscale = 1e-3

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  fa = ca+ma
  fb = cb+mb
  casa  = (/"A: "+na, "B: "+nb/)
  caseM  = (/pa+fa,pb+fb/)
  case  = (/caseM(0),caseM(1)/)
  Ncase = dimsizes(case)
  plotname = "figure_"+vname+"_260E-slice_"+shortname

  print(" ")
  print(" ")
  print("..................................................................")
  print(" ")
  print(" ")
  print(" varname   : "+plotname)
  print(" ctrl      : "+case(0))
  print(" test      : "+case(1))
  print(" plot name : "+plotname)
  print(" ")
  print(" ")
    
;************************************************
; Open netCDF file and read in dimensions
;************************************************

   filename=case(0)+".nc"
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)
   print("nlev "+sprinti("%7.2i",nlev))
   print("nlat "+sprinti("%7.2i",nlat))
   print("nlon "+sprinti("%7.2i",nlon))
   
;************************************************
; set default res
; res1 for gsn_csm_contour
; resp for gsn_panel
;************************************************
  wks   = gsn_open_wks (plotfmt, plotname)             ; open workstation
  gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")          ; choose colormap
  ;gsn_define_colormap(wks,"wh-bl-gr-ye-re")          ; choose colormap
  res1                      = True                ; plot mods desired
  res1@gsnFrame             = False               ; advance page
  res1@gsnDraw              = False                
  res1@gsnSpreadColors      = True                 ; use full colormap 
  ;res1@gsnSpreadColorStart  = -40
  ;res1@gsnSpreadColorEnd    = 40
  res1@cnFillOn             = True                 ; turn on color
  res1@cnLinesOn            = False                ; turn on/off contour lines
  res1@cnFillMode           = "AreaFill"           ; default "AreaFill" "RasterFill"
  res1@lbLabelBarOn          = True               ; turn on/off the labelbar
  res1@lbLabelsOn            = True                ; turn on/off the label of labelbar
  res1@lbOrientation         = "Vertical"          ; labelbar Orientation "Vertical"
  res1@lbBoxMinorExtentF     = 0.40                ; wide of the LabelBar 
  res1@lbLabelFontHeightF    = 0.016                                               
  res1@lbLabelAutoStride     = False                ; optimal labels         
  res1@cnLineLabelsOn        = False                ; turn on line labels
  res1@cnLinesOn             = False                ; turn on line 
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  res1@cnLevelFlags         = False
  res1@tiXAxisFontHeightF    = 0.020               ; XAxis Font
  res1@tiYAxisFontHeightF    = 0.018               ; YAxis Font
  res1@tmXBLabelFontHeightF  = 0.015               ; YXxis Font 
  res1@tmYLLabelFontHeightF  = 0.015               ; YAxis Font 
  res1@tiYAxisString         = "Pressure(hPa)"
  ;res1@tiXAxisString         = "Latitude"
  res1@vpXF                  = 0.12     ; x location left
  res1@vpYF                  = 0.88     ; y location  up
  res1@vpWidthF              = 0.55     ; width
  res1@vpHeightF             = 0.30     ; height
  res1@gsnStringFontHeightF  = 0.018
  res1@tiMainString          = " "     ; title    
  res1@tiMainFontHeightF     = 0.03
  res1@tiMainFontThicknessF  = 0.35
  res1@trYReverse            = True
  res1@gsnYAxisIrregular2Linear = True
  resP                 = True          
  resP@gsnPanelYWhiteSpacePercent = 1  ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 3  ; draw panel with white space added
  resP@txString            = ""
  resP@gsnPanelLabelBar    = False               ;  common colorbar
  resP@lbOrientation       = "Horizontal" 
  resP@lbLabelFontHeightF  = 0.016               ; labels size
  resP@lbBoxMinorExtentF   = 0.30                ; area allotted to each box of the LabelBar 
  resP@lbLabelAngleF       = 15
  res1@tmYRMode             = "Automatic"          ; turn off special labels on right axis

  
;************************************************
; read in variables & create plot
;************************************************

  plot = new(Ncase*3,graphic)

  do i=0,Ncase-1
     filename=case(i)+".nc"
     fi=addfile(filename,"r")
    ;print(" read "+case(i))

     freqi = fi->FREQI
     varx  = fi->IN_hom
     vary  = fi->IN_het 
     PS    = fi->PS

     freqi = where(freqi.lt.1.e-5,-999,freqi) 

     fvarx = varx 
     fvary = vary

     fvarx = (/ varx/freqi /) 
     fvary = (/ vary/freqi /) 

     printVarSummary(varx)  
     printVarSummary(fvarx)  
     
     P0mb =1000.
     hyam = fi->hyam
     hybm = fi->hybm

     ; type of interpolation: 1 = linear, 2 = log, 3 = loglog
     interp = 2
     ; is extrapolation desired if data is outside the range of PS
     extrap = False

     ; create an array of desired pressure levels:
    ;pnew = (/ 1000.,925.,850.,800.,750.,700.,600.,550.,500.,450.,400.,350.,300.,250.,200.,150.,100./)

     ;; use cam5 l30 levels, find and use the standard levels 
     
     pnew = (/ 100., 121., 143.,  \ 
               168., 200., 233., 274., \ 
               300., 379., 446., 500., \ 
               610., 700., 763., 821., \ 
               850., 887., 913., 936., \ 
               957., 976., 993., 1000. /) 
 

     pnew = pnew(::-1) 

     nlev2 = dimsizes(pnew) 
     
     ; note, the 7th argument is not used, and so is set to 1.
     pvarx = vinth2p(fvarx,hyam,hybm,pnew,PS,interp,P0mb,1,extrap)
     pvary = vinth2p(fvary,hyam,hybm,pnew,PS,interp,P0mb,1,extrap)

     printVarSummary(pvarx)  

    ;var3dx = dim_avg_Wrap(pvarx(time|:,lev_p|:, lat|:, lon|:)) 
    ;var3dy = dim_avg_Wrap(pvary(time|:,lev_p|:, lat|:, lon|:)) 
     var3dx = pvarx(:,:,:,{260.})
     var3dy = pvary(:,:,:,{260.}) 

     var3dx!0 = "time"
     var3dx!1 = "lev_p"
     var3dx&lev_p = pnew
     var3dx&lev_p@units = "hPa"
     var3dx!2 = "latitude"

     var3dy!0 = "time"
     var3dy!1 = "lev_p"
     var3dy&lev_p = pnew
     var3dy&lev_p@units = "hPa"
     var3dy!2 = "latitude"
     
     var2dx = var3dx(0,:,:) 
     var2dy = var3dy(0,:,:) 
     
     var2dx@_FillValue = -999
     var2dy@_FillValue = -999

    ;printVarSummary(var2d) 
    ;print(var2d&lev_p@units) 

     var2dx=var2dx*vscale
     var2dy=var2dy*vscale
     
     res1@gsnLeftString = casa(i)
     res1@gsnRightString= vunits                


    ;print(var2dx) 

     var2dx = where(var2dx.le.0.,-999,var2dx) 
     var2dy = where(var2dy.le.0.,-999,var2dy) 
 

     res1@gsnSpreadColors      = True                 ; use full colormap 
     res1@gsnSpreadColorStart  = 20                    ; start at color 10
     res1@gsnSpreadColorEnd    = 195                    ; end at color 96
  
     res1@gsnRightString= "#/L" 
     delete(res1@cnLevels)
     res1@cnLevels   = (/0.1,0.5,1,2,5,10,20,50,100,200,500,1000/)  
     
     ;;.............................................................................................
     ;; 1st panel 
     ;;.............................................................................................
     
     n=i+0
     res1@gsnLeftString = casa(i)
     res1@gsnCenterString = "HOM"
     plot(n) = gsn_csm_pres_hgt(wks,var2dx(:,:), res1) 

     ;;.............................................................................................
     ;; 2nd panel 
     ;;.............................................................................................
     
     n=i+2
     res1@gsnLeftString = casa(i)
     res1@gsnCenterString = "HET"
     plot(n) = gsn_csm_pres_hgt(wks,var2dy(:,:), res1)
     
     ;;.............................................................................................
     ;; 3rd panel 
     ;;.............................................................................................
     
     do nz=0,nlev2-1
     do ny=0,nlat-1
        if (.not.ismissing(var2dy(nz,ny)).and..not.ismissing(var2dx(nz,ny))) then 
           if ((var2dy(nz,ny)+var2dx(nz,ny)).ge.0.5) then  
              var2dy(nz,ny)=var2dy(nz,ny)/(var2dy(nz,ny)+var2dx(nz,ny)) 
           ;else
        ;  var2dy(nz,ny)=-999  
           end if
        else
           if(.not.ismissing(var2dy(nz,ny)).and.ismissing(var2dx(nz,ny))) then 
              var2dy(nz,ny)= 1. 
           else
              var2dy(nz,ny)=-999 
           end if 
        end if
     end do
     end do
     
     delete(res1@cnLevels)
     res1@cnLevels   = (/1,10,20,30,40,50,60,70,80,90/) 
     n=i+4
     res1@gsnLeftString = casa(i)
     res1@gsnCenterString = "               HET/(HOM+HET)"
     res1@gsnRightString= "%"  
     var2dy=var2dy*100
     plot(n) = gsn_csm_pres_hgt(wks,var2dy(:,:), res1) 
     
     delete(var2dy)
     delete(var2dx)
     delete(var3dx)  
     delete(var3dy) 
     delete(hyam) 
     delete(hybm) 
     delete(freqi) 
     delete(varx) 
     delete(vary) 
     
     delete(pnew) 
     
  end do 
  
  ;;gsn_panel(wks,plot,(/3,4/),resP)
  gsn_panel(wks,plot,(/3,2/),resP)

  
end


